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Abstract 

This paper extends the gas-kinetic scheme for one-dimensional inviscid shallow water equations (J. 
Comput. Phys. 178 (2002), pp. 533-562) to multidimensional gas dynamic equations under gravi- 
tational fields. Four important issues in the construction of a well-balanced scheme for gas dynamic 
equations are addressed. First, the inclusion of the gravitational source term into the flux function is 
necessary. Second, to achieve second-order accuracy of a well-balanced scheme, the Chapman-Enskog 
expansion of the Boltzmann equation with the inclusion of the external force term is used. Third, 
to avoid artificial heating in an isolated system under a gravitational field, the source term treat- 
ment inside each cell has to be evaluated consistently with the flux evaluation at the cell interface. 
Fourth, the multidimensional approach with the inclusion of tangential gradients in two-dimensional 
and three-dimensional cases becomes important in order to maintain the accuracy of the scheme. Many 
numerical examples are used to validate the above issues, which include the comparison between the 
solutions from the current scheme and the Strang splitting method. The methodology developed in 
this paper can also be applied to other systems, such as semi-conductor device simulations under 
electric fields. 

1 Introduction 

Since most astrophysical problems are related to the hydrodynamical evolution in a gravitational 
field, the correct implementation of the gravitational force in an astrophysical hydrodynamical code is 
essential, especially in systems with long time integrations, such as modeling star and galaxy formation. 
Even though many hydrodynamical codes have been successfully applied to astrophysical problems, 
including the Piecewise parabolic method (PPM) and total variation diminishing (TVD) codes [5, 
12], most have considered only short time evolutions with strong shock or expansion waves. With 
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the slowness of galaxy evolution, many codes have difficulties due to the improper treatment of the 
gravitational force effect, the so-called source term in the Euler or the Navier-Stokes equations. After 
integration of millions of time steps, many codes cannot settle down to determine an isothermal 
steady state solution for a gas in a time-independent external gravitational field. The solution will 
either oscillate around the equilibrium solution, or simply deviate from equilibrium due to artificial 
heating, which triggers numerical gravitation-thermal instability, i.e., the collapse of the gas core. 
There have been many attempts to construct such a well-balanced gas dynamic code that preserves 
the hydrostatic solution accurately [8, 20, 2]. 

In recent years, the study of flow equations with source terms has attracted much attention in the 
computational fluid dynamics community [18]. Flow equations with source terms in one-dimensional 
(ID) case can be written as 

(1) Wt + F{W), = S, 

where W is the vector of conservative flow variables with the corresponding fluxes F{W). Here, S is 
the source, such as the gravitational force. If we integrate the above equation with respect to dx in a 
spatial element, j, from Xj_i/2 to Xj_^_l/2, and dt in a time interval from to t'^'^^, 

/ / {Wt + F{W),;)dxdt= / / S{W)dxdt, 

we get 

3 1/2 

where Ax = 2:^+1/2 ~ is the cell size and Wj is the averaged value of W in cell j, i.e., 

T^, = — / Wdx. 

^X Jxj_i/2 

The above equation is an integral form of Eq.(l) and is exact and equivalent to Eq.(l). In an earlier 
paper on the shallow water equations [16], we emphasized the importance of including the source term 
effect into the flux function at a cell interface in order to develop a well-balanced scheme. However, 
the exact equilibrium solution, such as u = and h + B = constant, for the shallow water equations, is 
relatively simple and special. The simplicity can be realized through the relation dh/dx = —dB/dx = 
constant inside each computational cell once the bottom topology is approximated as a linear function 
inside each cell. Therefore, the MUSCL-type scheme with a single slope inside each cell can be precisely 
used to reconstruct the initial well-balanced solution. However, these techniques developed for the 
shallow water equations cannot be directly extended to gas dynamics equations. For example, the 
hydrostatic equilibrium solution for the gas dynamic equations may become p ~ e~"^, which is an 
exponential function. This exact equilibrium solution cannot even be precisely reconstructed as an 
initial condition in the MUSCL-type scheme with a single limited slope inside each numerical cell. 
Even though, it seems impossible to develop a well-balanced scheme for the gas dynamic equations 
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with an error up to the order of machine zero, a proper numerical treatment is still essential to get 
accurate solutions. 

For a gas flow under an external time-independent gravitational field, there exists an isothermal 
solution, the so-called hydrostatic or well-balanced equilibrium solution, with a constant temperature 
and zero fluid velocity. This solution is an intrinsic solution due to the balance between the flux 
gradient and source term, i.e., F(W)x = S. However, based on the macroscopic equations, it is 
unclear how the gas settles down to an isothermal solution. Usually, additional assumptions, such 
as a constant temperature, are needed before deriving such an equilibrium result. Numerically, it is 
challenging to obtain such a solution. If Eq.(l) is solved by simple operator splitting method, such as 
Wt + Fx = for the flux at a cell interface and Wt = S to account for the source term inside each 
cell, as observed frequently, many schemes simply deviate from the exact solution after millions of 
time step integrations. The systematic errors will accumulate and eventually crash the calculation. 
This happens to be the so-called numerical heating phenomenon, which was reported in [13]. Due 
to the accumulation of numerical errors, the internal energy in an isolated gas system with adiabatic 
boundary conditions under a constant gravitational field will increase forever. Figure 1 shows such 
a solution, where the internal energy goes up after millions of time steps. The non-conservation of 
the total energy in the system is not surprising because the gas dynamic equation (1) itself is not 
written in a conservative form. The problem in the gas dynamic equations is that the numerical error 
accumulates in one direction. On the contrary, there is no such problem in the shallow water equations 
even if the numerical scheme is not a well-balanced one. 

Based on the Boltzmann equation, a flow under a gravitational field in a hydrostatic state can be 
described as 

^.^_^.^ = 

dx dx dc ' 

where / is the gas distribution function, c is the particle velocity, and (f) is the external gravitational 
potential. The general solution of the above differential equation is 

/(f,c)=.^(0+^c2), 

where T is an arbitrary function of = c - c = + + c^. This equilibrium solution clearly 
shows that the gas will have the same temperature, because the temperature becomes a constant 
multiplier for the function (j) + ^c?. For example, the function J- in hydrostatic equilibrium will 
become exp(— (^c^ -|- (f))/kT), where kT is a constant temperature. The above analysis shows the 
usefulness of using a kinetic equation instead of macroscopic equations to describe a gas flow under a 
gravitational field. It also implicitly demonstrates that the flux functions, which are the moments of 
/, have to take into account the source term effect, i.e., (p, in the solution of /. 

In the past several years, a gas-kinetic BGK Navier-Stokes (BGK-NS) scheme has been developed 
for compressible Navier-Stokes equations [15, 9, 10]. It achieves success in computing viscous flow 
solutions, such as the hypersonic viscous and heat conducting flows [14]. Due to the simplicity of 
particle transport, a multi-dimensional scheme can be constructed by incorporating a multidimen- 
sional particle propagation mechanism through the cell interface, where the particle transport in both 
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normal and tangential directions of a cell interface can be taken into account. However, constructing 
a well-balanced BGK-NS method for gas dynamic equations under a gravitational field is not straight- 
forward. As shown above, with a hydrostatic state, the density inside each numerical cell is distributed 
exponentially in space, which cannot be precisely reconstructed by a simple limited slope in a standard 
MUSCL-type scheme. In other words, even starting from the initial condition, some numerical error 
is already introduced. Another important issue that needs to be addressed is the Chapman-Enskog 
expansion of the Boltzmann equation, where the external force effect has to be added to account for 
the balance between the dissipative and external forcing terms. This is important to capturing the 
viscous flow solution, especially with a scheme that should work under both continuous and discontin- 
uous flow conditions. To achieve a well-balanced scheme, the gravitational force effect on the particle 
acceleration and on the numerical flux at the cell interface needs to be explicitly taken into account 
as well. In order to avoid the numerical heating (cooling) effect, a consistent treatment between the 
source term inside each cell and the flux at the cell interface is required. The BGK scheme presented 
in this paper is a well-balanced scheme up to second-order accuracy. On the contrary to the shallow 
water solutions, the hydrostatic equilibrium state for the gas dynamic equations, such as zero flow 
velocity, cannot be precisely kept by the current scheme up to the machine zero, i.e., 10~^^. How 
to extend a well-balanced strategy from the shallow water equations to the gas dynamic equations 
remains an open problem. 

This paper is organized as follows. Section 2 gives the details of the multi-dimensional BGK-NS 
scheme under a gravitational field. Section 3 presents the results from many test cases, which confirm 
the importance of some strategies in the development of a well-balanced scheme. At the same time, 
the numerical heating phenomenon is presented. The functions from the nonlinear limiter and the 
multi-dimensional transport mechanism on the well-balanced scheme are described. The last section 
is the conclusion. 

2 BGK scheme in three-dimensional space 

In this section, the multidimensional finite volume gas-kinetic scheme based on the Bhatnagar-Groos- 
Krook (BGK) model for the compressible Navier-Stokes equations [15] is extended to solve a flow 
system under an external gravitational force. The governing equations solved by the gas-kinetic 
scheme are written as, 

(2) I- / Udn+ </F-dS= / QdO, 

ot Jn Js Jn 

where is the control volume, S is the surface of this volume with its surface normal vector pointing 
outward. The macroscopic variables are 

U{x,t) = ip,pu,Ef, 
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where p is the density, u is the velocity, and E is the sum of the thermal and kinetic energy. The 
fluxes take the form. 



F{x,t) 



^ pu ^ 

pu^u + pi — Tj 
^ {E + p)u - kVT -t-u j 



where (8) denotes the tensor product of the vectors, k is the thermal conductivity coefficient, I is the 
unit tensor, and E stands for the viscous shear stress tensor, 

= p{diUj + djUi) - ?(V • u)5ij, 

where /x and are the first and second viscosity coefficients. The source terms can be written as 
follows: 

Q = iO,-pV4>,-pu-V4>f, 

where (j) = 4>{x, t) is the gravitational potential. 

The BGK scheme is a finite volume method. In Cartesian coordinates, equation (2) is discretized 

as ^ ^ 

(3) C7;+1 = C7; - ^— /■* ^ (i? +v2(t) - F,-i/2(i))dt + r ^ Qjdt, 

where Uj and Qj represent the cell-averaged quantities in the jth cell, and -Pj+i/2(i) is the vector of 
time-dependent fluxes across the cell interface between the jth and j+lth cells. 

Generally, there are three stages in a shock capturing finite volume scheme, i.e., reconstruction of 
initial data, gas evolution, and projection. In the reconstruction stage, a piecewise continuous flow 
distribution inside each cell is obtained. Then, in the gas evolution stage, a time-dependent flux, 
Fj^i/2{t), is computed based on the local solution at the cell interface, such as the solution of a 
Riemann problem. Based on the fluxes, the cell average flow quantities can be updated to a new time 
level. At the same time, the source term contribution has to be added on the time evolution of the 
flow variables inside each cell. 

2.1 Reconstruction 

In the higher order finite volume method, the cell-averaged flow variables have to be reconstructed 
first to obtain the quantities at the cell boundaries. The reconstruction depends solely on the flow 
distribution at the beginning of each time step. In the smooth flow region, a direct connection of cell 
averages across a cell interface is a good choice. 



(4) C/,-+i/2 



Uj+i + Uj 



where Uj+1/2 stands for the flow variables at the interface, and the slope across the cell interface 
becomes {Uj+i — Uj)/{xj-\-i — Xj). A higher-order interpolation can be used to improve the spatial 
accuracy of the scheme. In a discontinuous flow region, the kinematic dissipation is added implicitly 
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to the reconstructed initial data by converting kinetic energy into thermal one [15]. The amount 
of kinematic dissipation depends on the limiters used. In the current study, the van Leer limiter is 
adopted. For the sake of simplicity, let Uj+i{xj^i/2) denote the reconstructed data on the right-hand 
side of the interface {xj^i/2) and let Uj{xj_^_i/2) be the data on the left-hand side. Then, 

(5) Uj{xj+i/2) = Uj{xj) + (xj-+i/2 - Xj)Lj, 

(6) C/j+i(xj+i/2) = C/j+i(a;j+i) - (x^+i - Xj+i/2)Lj+i. 
In the above expressions, the van Leer limiter is defined as, 

where s_) = sign(s_|_)-|-sign(s_), s_|_ = {Ujj^i — Uj)/{xjj^i—Xj), and s_ = {Uj — Uj-i)/{xj—Xj^i). 

2.2 Gas evolution 

In the BGK scheme, the flow evolution is governed by the BGK equation [1], 

(7) ^ + c- • V/ + (-V0) ■ = 

where f{x,c,t) is the gas distribution function, c is the particle velocity, r is the collision time, and 
= {d/dci,d/dc2,d/dc3). The right-hand side of (7) is the so-called relaxation model, which is 
a simplification of the collision term of the Botlzmann equation [19]. The equilibrium state, g, in 
equation (7) is a Maxwellian distribution, 



where ^ has N internal degrees of freedom, such as the molecular rotation, A is related to the gas 
temperature by A = m/2kT, m is the molecular mass, k is the Boltzmann constant, and T is the 
temperature. The above expression is for a gas in three dimensions. For ID and 2D flow simulations, 
the degree of freedom K of the variable is defined as K = N + 2 for ID flow and K = N + 1 for 2D 
flow, where 2 stands for the particle random motion in the X2 and X3 directions and 1 accounts for 
the X3 direction random motion only. In a 3D simulation, K is equal to N. In the equilibrium state, 
the internal energy variable, is defined as = + .^l + ■ • • + ^k- 

The relation between the macroscopic variable, U, and the distribution function, /, is 

"+00 ^ ^ ]_ 



r+oo ^ ^1 

U= V/dS, V=[l,c,-(c-c + aF, 
J —00 ^ 



where dS = dcd^ is the volume clement in the phase space with d^ = d4id42 • • • d^K- The moments of 
a Maxwellian are given in Appendix A of [15]. 

The Navier-Stokes equations can be regarded as a low-order approximation of the BGK equation 
for the linearly distributed flow variables in space with non-zero velocity and temperature gradients. 
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It has been shown that for a monoatomic [19] and polyatomic gas [17], the Navier-Stokes equations 
can be derived from the BGK model with the dissipative coefficients, 

2 AT N + 5 k 

3 A* + 3 2 m 

which are the dynamical viscosity coefficient, the second viscosity coefficient, and the thermal conduc- 
tivity, respectively. For an ideal gas, the equation of state becomes p = p/2X. 

The right-hand side of equation (7) describes a relaxation process from an initial non-equilibrium 
state, /, to an equilibrium state, g, through particle collisions. Here, r is called the particle relaxation 
time. Assume that r is a local constant and integrate equation (7) along its particle trajectory, 

d^' . ^ , , dc' ^ , 

We have 

(8) f{x, c,t) = - f 5e-(*-*')/Mt' + e-*/-f(xo, cq, 0), 

where < t' < t and 

co = c-(-V0)t, xo = x-ct+ -{-V4>)t^ 

are respectively the particle velocity and trajectory with initial position xq and velocity cq. The 
gravitational effect has been included in the particle acceleration and trajectory change. Since the 
term ^(— V^)t^ has a third-order effect on a numerical scheme, it will be ignored in the construction 
of the time-dependent solution, /, in the following. 

The integral of the right-hand side of equation (8) is generally difficult to evaluate because the 
equilibrium state, g, itself is implicitly related to /, and the initial distribution, f{xo,co,t' = 0), is 
unknown. In the BGK scheme, both are modeled as first-order expansions of a local equilibrium 
distribution with the inclusion of the Chapman-Enskog expansion for the dissipative term. To include 
the effect of gravitational force in the non-equilibrium gas distribution, the initial distribution, /o = 
f{xo,co,t' = 0), is further developed as 



(9) /o 
where 



gl[l -a!- -ct-V- ■ (-V0)t - r(a' -c + d ■ (-V^) + A% x < 
gl[l-ar ■ct-h' ■ {-V(t))t - TiaT ■ c + 6^ • (-V(/)) + A^)], x > 0, 



iV+3 



and 



N+a 

2 



are the equilibrium distributions at the left- and right-hand sides of the cell interface. The slopes in 
equation (9) are defined as a' = Vln^Q, (T = Vln^Q, 1? = Vglng^, ¥ = Va'lnt/Q, = dlngQ/dt, 
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= d In gydt. The above model is an extension of the Chapman-Enskog expansion in [15], where 
the gravitational force term has been added in the modification of the non-equilibrium state. The 
term proportional to r corresponds to the dissipative terms in the Navier-Stokes equations. The 
gravitational force makes two contributions on /q. One is the particle acceleration in the transport 
process to a cell interface, and the other is in the non-equilibrium initial state. Note that in the above 
initial gas distribution function, the flow variation in the tangential direction is also included in the 
definition of a} and (F. 

All slopes in (9) can be uniquely evaluated from the reconstructed flow variables and their deriva- 
tives. In the normal direction of a cell interface, {xj_^.i/2), the slopes on the left-hand side can be 
expressed as 

a\=^ ■ if, 

where a' is a local constant vector. Based on the relation between the distribution function and the 
macroscopic variables, we have 



Uj{Xj+i/2)-Uj{Xj) 
Xj+l/2 - Xj 



j a\^gi6E = j if V^g^dS. 



The above linear system can be solved directly for Qj and the detailed formulation is given in the 
appendix for the 3D case. Similarly, the slopes on the right-hand side of the interface can be obtained 
by solving 

Uj+ijxj+i) - ?7j+i(xj+i/2) _ 

Xj+l — Xj^i/2 

In the tangential direction (xk) of the cell interface, the slopes are obtained from 



Uj,k+l{Xj+l/2) - Uj,k-l{Xj+l/2) ^ 
Xj+l/2,k+l ~ Xj+l/2,k-l 

C^j+i,fc+i(a;j+i/2) - Uj+i^k-i{xj+i/2) 



Xj+l/2,k+l — Xj+i/2,k-l 

The terms multiplied by r in equation (9) represent the non-equilibrium part of a Chapman-Enskog 
expansion, and make no direct contribution to the conservative flow variables. Therefore, based on 

jict.c + P- {-V<f>) + A'WodE = 0, 
J (5- . C-+ . (- V0) + A^),Pgl,dE = 0, 

the temporal slopes A' and A^ can be determined. Note that, for the first time, the gravitational force 
effect on the determination of the non-equilibrium state has been proposed in the above equations. 

The equilibrium state g{x' ,(f ,t') in the integral solution is constructed by a Taylor series expansion 
of a local equilibrium distribution around {x,c,0). With the inclusion of the external source effects. 
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it has the form 



(10) 9 = 



go[l - (t - t')aLcn - (i - t'fbU-Vcj))n 

-{t- t'fat -ct-it- t')ht ■ (-V(/.)f + At'], c„ < 

50 [1 -{t- t')(f^Cn -{t- t')X{-V(f>)n 

-{t- t')~dt -ct-it- t')ht ■ {-Vct))t + At'], c„ > 



where 

50 = Po ( 

The parameters in (10) with subscript n represent the normal components of the vectors to the cell 
interface, and those with subscript t stand for the tangential components. This is again an extension 
of g in [15]. It will be shown in Section 3 that the inclusion of both the external source term and the 
flow gradients in the tangential direction in the expressions of /o (9) and g (10) is very important to 
obtain a well-balanced scheme and to improve the accuracy and stabihty of the scheme. 

In order to determine all parameters in g (10), the equilibrium state {Uq) at the cell interface at 
t = must be obtained first. Since the /o is totally determined, taking the limit t in equation 
(8) and substituting its solution into conservation constraint J{g — f)ipdE = at the cell interface, we 
have 



/ 



go^fdE = Uo = / / g^o^fdE + f f g^V^dS. 

Jc„>OJ Jc„<OJ 



ICn>OJ JCn<0. 

Then, the slopes in equation (10) can be obtained through 

Uo - Uj{xj) 



in the normal direction on the left- and right-hand sides of the cell interface, and 



/ attfdE = / algoV'dS + / a^goV^dS, 

J JCn>0 J JCn<0 J 



in the tangential direction. 

The only unknown in the expression of g (10) is A. By substituting the approximations (9) and 

(10) into the general solution (8), we have 

/ = {(1 - e-*/-) + (e-*/-(t + r) - T)[(aic„ + 6i(- V</.)„)F[c„] 

+(a;c„ + 6;(-V0).„)(l - H[c,,]) + cl^fct + bf (-V(/.)t] + [t - r(l - e-'/^)]A}go 
+ e-*/^{[l -af-ct-P- {-V(l))t - T(a} -c + P- (-V(/)) + A^)]H[cn]gl) 

(11) + [l-(f -ct-b' ■{-Vcf>)t-T{(f -c + F ■{-Vct>)+A')]{l-H[cn])g'o}, 



9 



where the Heaviside function H[x] is defined by 

H[x] = I 



0, a; < 0, 

1, x > 0. 



Due to the conservation of the conservative flow variables during particle collisions, the moment of the 
right-hand side of equation (7) should vanish. At the cell interface, the compatibility condition gives 

rAt 

' {g - /)V'dtdS = 0, 



from which we get 

70 j MgodE 



where 



70 
71 

72 

73 
74 

75 



+(a;c„ + 6;(-V<^)n)(l - H[cn]) + ~dfct + hf {-V(t>)t) 

+ 73(i?[Cnbo + (l-^M)55) 

+ (74 + 75) (i?[Cn] (3' • C + ? • (-V0))5^ 

+{l-H[cn]W-B+h^-{-Vc^))9l) 
+ -i^{H[cn]A'g'^ + (1 - H[cn])A'9l)]^d'E, 

= Jr{l- e-*/ndt = r[At - t(1 - e"^*/^)], 
-*Adt = -r[e-^*/"-l], 



= J -(e-*/^(t + r) - r)dt = - At - (74 + 75), 



-71, 



J te-*/^dt = -r[Ate-^*/^ - r(l - e'^^^^)], 
j re~*/'^dt = r7i. 



Therefore, A is fully determined from the above equations. As a result, the time-dependent gas 
distribution function, f{x,c,t), is obtained. The corresponding fluxes in the Xj direction across the 
cell interface can be computed by 



= j [^i9o + To^^o 

+ T2i{aicn + bU-Vcl>)n)H[cn] 

+(a;c„ + 6;(-V0)n)(l - H[cn]) +SfCt + %- {-V(f))t)go 

+ Ts{H[Cnyo + il-H[CnM) 
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+ (r4 + T^){H[cn]{n' -c + ff- {-Vcj)))gl, 

+(l-i7[c„])(cf •c + 6''.(-V0))5S) 
(12) + T,iH[cnWgl + (1 - H[cn])A^g'oMcjdE, 



where 



ro = 


2 


Fi = 


At + T 


r2 = 


-T72, 


rs = 


-T73, 


r4 = 


-■7"74, 


rs = 


-T75- 



The above formulation is based on the discontinuous initial condition at the cell interface, and the 
above scheme can be faithfully applied to a flow with discontinuities. For a smooth flow, based on the 
reconstruction without using the limiter, the cell interface discontinuity will automatically disappear. 
So, under this condition, we have 

di = gr = S, 9 = b'' = b, A^ = A^ = A. 

Therefore, the solution of the BGK equation (8) will simply go to 

(13) fix,c,t) = g{x,c,0)[l-T{S-c + b-{-V(f>) + A)+At]. 

Prom this solution, we can clearly see that the BGK solution from a general discontinuous initial 
condition will go to the well-defined continuous solution automatically, i.e., the Chapman-Enskog 
expansion for a smooth flow under a gravitational field. In other words, the multidimensional BGK 
scheme satisfies the consistency requirement. Within a time scale, r, the particle acceleration due to 
the gravity will distort the distribution from the equilibrium through the term b ■ (— V(/)). 

It is well known that the BGK model has an intrinsic unit Prandtl number. To make the BGK 
scheme simulate any heat-conducting flow, the energy flux related to the heat conduction part in (12) 
can be modified [15] to 

rr"" = Fe + (1/Pr - 1)^, 

where q is the heat flux in the normal direction of the cell interface and Pr is the Prandtl number. 
According to the gas-kinetic theory, the heat flux is a result of energy transport carried by molecules 
in their random movement. In the Xj direction, q can be evaluated by 

Qi = ^Jici-Ui){{c-u)-{c-u)+^'^)fdE. 

With a discontinuous initial condition, the collision time has the form [15], 

PO \Pl/k + Pr/K\ 
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where the first term comes from the Chapman-Enskog theory of the BGK model, which states that 
the particle collision time is related to the ratio of the dynamical viscosity coefficient to the pressure. 
The second term enhances the dynamical dissipation in cases of flow discontinuity. In the smooth flow 
region or around the slip line, the numerical dissipation introduced above is very small because of the 
continuous pressure distribution. 

The above formulation clearly shows the complexity of the gas- kinetic scheme for the Navier-Stokes 
solutions under a gravitational field. The splitting of a Maxwellian distribution into two half ones leads 
to a very complicated flux splitting method with error functions and terms involving exponentials. The 
drawback may be overcome by Pcrthamc's compactly supported equilibrium distribution function [11] 
instead of a full Maxwellian distribution. This will makes BGK scheme much simpler and attractive, 
which deserve further investigation. 

2.3 Projection 

In the projection stage, the average mass, momentum, and energy inside each cell are updated. Any 
inconsistent numerical treatment between the flux and the source term in this step may easily cause 
numerical heating in the simulation of gravitational gas dynamics. Here, we are going to analyze 
the so-called numerical heating reported by Slyz [13]. When numerical heating occurs, the computed 
system will continuously increase its internal energy, which is shown in Fig. 1. 

Typically, there are two ways to update the energy inside each cell. One is the non-conservative 
form: 

at Ja Js 

(14) = I pu- {-V4>)dn. 

For a time-independent potential, the gravitational energy can be absorbed into the total energy. 
Therefore, the energy equation in this special case can be written in a conservative form as the 
following: 

(15) ^ / Etotdn + / ((Etot + P)n - kVT - S • u) • dS = 0, 

ot Jn Js 

where Etot = E + pcf). The discretized forms of equations (14) and (15) are 

, „fn+l ^n + l 

(16) = / {Fe{xj+,/2, t) - Fe{xj_,/2, mt + / QE(xj)dt, 

(17) EZt\xj) = Eltixj) - ^— /* {FE,^,{Xj+y2,t) - FE,A^j-i/2,t))dt, 

^j+l/2 - Xj-l/2 Jf^ 

respectively, where the gravitational source term is 



J QE{xj)dt = ^J J^pu- {-V(p)dndt. 



Slyz and Prendergast concluded that equation (17) rather than (16) should be used to update the 
energy term. Otherwise, numerical heating will be introduced into the numerical solution because of 
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the inaccurate computation of the source term [13]. Eventuahy, the numerical gravitational system 
will break down. This argument is perfectly correct. However, it is not totally complete. With a 
careful discretization of the source term in (16), we can still establish a stable and accurate solution. 
Numerically, the inaccuracy comes mainly from the low-order approximation and the inconsistent 
computation of the source term. Here, inconsistency means that the same flow variable at different 
locations of the governing equations are calculated by different methods. For example, in the pro- 
jection formulation, if the gravity is time-independent, then the source term in equation (16) can be 
approximated by 

(18) flJJn'^' (-'^'^)'i^dt = ^J J pudt • (- V0)dJ7, 

where the term in the square brackets is the time integral of the mass flux, which should be the same 
as the corresponding term in the continuum equation. In Slyz and Prendergast's paper [13], the mass 
flux in the continuous equation was calculated by an early BGK scheme (denoted as F^^^), while the 
mass flux in equation (18) was computed by a bilinear interpolation (denoted as i^*). Suppose their 
difference is — F^^^ = SQ. If all the fluxes in (3) are calculated by the BGK method while the 
source term in (16) is calculated by F^* • ^ / {—'V(f))dQ,, the additional source term, (^Q • ^ /(- V^^)dO, 
will appear in the energy equation, (16). This could be the origin of numerical heating. Unfortunately, 
the error will grow in one direction in a system under gravity. The instability is closely related to the 
so-called gravo-thermal instability. But the reason is due to the numerical conservation error, instead 
of the physical one due to the energy loss in an astrophysical system. 

For an isolated or adiabatic system, after a long time integration, any inconsistent treatment of the 
source term in the formulation (16) will generate numerical heating and the system will deviate from 
its physical solution. In order to avoid the instability and inaccuracy in the solution, we have taken 
the following two steps when using (16). The first one is to compute the source term consistently. In 
the current study, the mass flux at the cell interface is used to approximate the integration inside the 
cell. The source term, (18), is discretized as 

J QE{xj)dt = \j (Fjgk(xj_,/2) +F5g'^(xj+i/2))(-V^)jdt. 

The second one is that the gravitational forcing effect is explicitly included in the flux evaluation 
as presented for the BGK scheme in the last subsection. Equipped with these two mechanisms, the 
scheme with (16) is still a well-balanced method up to the second-order accuracy. 

When equation (17) is used, the time-integrated total energy flux can be written as 

Fe^., =Fe + Fp4>. 

Therefore, the total energy will be precisely conserved. Even through the gravitational force effect is 
not included in the flux evaluation in [13], based on the formulation (17) the numerical error will not 
grow and accumulate forever. The departure from the isothermal hydrostatic state will stay bounded 
due to the total energy conservation, but with large oscillations. Even when using (17), the numerical 
accuracy can be much improved once the gravitational force is included in the flux construction. The 
numerical results presented in the next section support the above argument. 
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3 Numerical tests and discussion 



In this section, numerical tests of ID, 2D and 3D cases are performed to validate the gas-kinetic 
method. Each of them is very sensitive to the accuracy of the scheme. Some of them are run for 
millions of numerical steps whereby the accumulation of any small error would become significant for 
such a long time integration. 

3.1 Perturbation of the ID isothermal equiUbrium solution 

This first test case is from the the paper by LcVcquc and Bale [8], where a quasi-steady wave- 
propagation algorithm was developed for an ideal gas subject with a static gravitational field. Initially, 
an ideal gas with 7 = 1.4 stays at an isothermal hydrostatic state, 

po{x) = poix) = e""^, and uq{x) = 0, 

for X G [0,1]. The gravity acts in the negative x direction, such as — V^(x) = —1. Initially, the 
pressure is perturbed by 

p(a;,i = 0)=po(a;)+r/e-"(^-^°)', 

where a = 100, xq = 0.5, and 77 is the amplitude of the perturbation. 

The computation is conducted with 100 grid points in the whole domain. The results from the BGK 
scheme with and without the gravitational source term included in the fluxes are shown in Figure 2. 
The figure clearly shows that the inclusion of the source term in the flux function is very important 
for the accurate capturing of small perturbations. In the above calculations, the BGK scheme uses 
the van Leer limiter in the initial reconstruction. Therefore, the BGK scheme itself works in both the 
smooth and discontinuous flow regions. However, the quasi-steady wave propagation algorithm in [8] 
works only in the smooth flow region. 

3.2 One- dimensional gas falling into a fixed external potential 

This case is taken from Slyz and Prendergast's paper [13] to investigate the numerical accuracy of the 
BGK scheme. The gas is initially stationary {u = 0) and homogenous {p = 1, e = 1, where e is the 
internal energy). The gravitational potential has the form of a sine wave, 

C . 2ttx 



where £ = 64 is the length of the computational domain and = 0.02. The ratio of the specific 
heat, 7, is 5/3. The periodic boundary conditions arc implemented in this system. After the initial 
transition, the system eventually reaches an isothermal hydrostatic distribution, where the tempera- 
ture settles to a constant (T{x,t) = Tq) and u = 0. For an ideal gas (p = pRT), from the hydrostatic 
equation dp/dx = p{—d(j)/dx), the density distribution becomes p{x,t) ~ exp {—(p/RTo). However, 
the system is highly nonlinear since Tq depends on p. It is thus better to use a numerical method to 
solve this problem. Each simulation result presented here is obtained with 500, 000 time steps on a 
Cartesian grid with Ax = 1. The Courant number used is 0.9. 
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When the numerical heating occurs in an inconsistent scheme, the gas will become hotter due 
to the numerical error. As used by Slyz, the schemes solving equation (14) are named the energy 
source term (EST) schemes, and others solving equation (15) are called the energy conservation 
(ECT) schemes. The final state (density, velocity and temperature) of the system given in Figure 3 
is calculated from both the EST and ECT BGK schemes. A consistent source term treatment for 
the EST scheme is used, such that the mass flux integration inside each cell is approximated by the 
average of the cell interface mass fluxes. In the above calculation, the gravitational effect is included in 
the flux evaluations for both schemes. In this figure, the diamonds indicate the solutions of the ECT 
scheme and the solid lines indicate the solutions of the EST scheme. Figure 3 shows that the results 
obtained by the EST scheme and ECT scheme arc almost indistinguishable. The small derivation 
from the hydrostatic state in both schemes is due to the inaccurate representation of the exponential 
function by a single slope inside each cell. We cannot find any simple way to remedy this. 

Even within the framework of the ECT scheme, in the following we are going to show that it is 
necessary to include the gravitational force in the flux for maintaining high accuracy. For a stationary 
state, n = 0, it can be shown that in the current scheme the fluxes induced by the spatial derivatives 
J a ■ ccigipfd'E are totally balanced by the gravitational force terms / b ■ {—'V(p)cig'ipfd'E, where there 
is no adjustable parameter. In other words, without considering the external force effect in the flux, 
the steady state cannot be accurately reached by a numerical scheme. The solution given in Figure 4 
is calculated by ECT scheme. The diamonds are the solution of the ECT scheme with a gravitational 
effect in the flux and the solid lines are the ones without the gravitational effect in the flux. Figure 4 
shows that neglect of the gravitational effect in the gas evolution stage will drive the solution away 
from the hydrostatic equilibrium solution even though the solution will not blow up due to the total 
energy conservation, the so-called ECT scheme. We believe that any scheme based on the Riemann 
solution, if the gravitational term is not included in the flow evolution, may have a similar problem. 
By including gravity, the relative fluctuation in both the velocity and temperature is much reduced. 
It is difficult to predict if there will be any shock capturing well-balanced scheme, where the fluid 
velocity in a hydrostatic equilibrium state can be kept up to the machine zero. 

This test case illustrates that a simple operator splitting approach does not work and the inclusion 
of the gravitational force in the numerical flux is the key for a well balanced scheme. However, it 
is widely recognized and proved mathematically that the Strang splitting method is a second order 
accurate time-integration method for general evolutionary equations with multiple operators. In the 
following, we apply the Strang splitting method to this test case as well. The update of the flow 
variables are 

(19) C/"+^ = S[At/2]HBGK[At]S[At/2]U'', 

where S represents the source term contribution inside each cell, and Hbgk is the BGK flux across a 
cell interface without including the gravitational force term. The results from the EST BGK scheme 
and the above EST Strang splitting method are shown in Figure 5. The Strang splitting method gives 
accurate velocity solution, which is similar to that from EST BGK. However, the temperature of the 
EST Strang solution is continuously increasing with the total number of integration steps. For a slow 
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gas evolution system under gravitational field, the numerical heating from the Strang splitting method 
will become severe. 

3.3 Onset of two-dimensional compressible convection 

Consider a fluid confined by two horizontal parallel plates, where the gravity is imposed vertically, and 
the upper plate is fixed at a lower temperature while the lower plate is fixed at a higher temperature. 
For inviscid fluids, when the temperature gradient exceeds the adiabatic temperature gradient, the 
system becomes convectively unstable. In the astrophysics literature, such a condition is known as 
the Schwarzschild criterion for convection ([7],p39). For general cases where the viscosity is taken into 
account, the motion is always impeded by the viscosity, and the temperature gradient is reduced by 
the thermal conductivity. To study the general situation, the Raylcigh number is used to indicate the 
combined effects of the dynamical viscosity, thermal diffusion and temperature gradient. The steady 
convection will not occur unless the Rayleigh number is larger than a specific value, i.e., the critical 
Rayleigh number. The compressible convection is a sensitive test on the accuracy and stability of 
any numerical scheme. The onset of stable large-scale circulation depends sensitively on the numerical 
dissipation and boundary treatment. A low-order scheme may introduce too much artificial dissipation, 
which increases the critical Rayleigh number. Due to the large scale in height and pressure changes, 
the instability can be triggered easily. 

In the gas-kinetic BGK scheme, the dissipative coefficient is controlled by the particle collision 
time. In a continuous flow region, only physical viscosity needs to be used. However, there are two 
ways to add the artificial dissipation in a numerical scheme to enhance the shock-capturing ability. The 
kinematic dissipation can be introduced through limited reconstructed initial data, which implicitly 
transfers kinetic energy into thermal energy. The dynamical dissipation can be added as well by 
enlarging the viscosity coefficient in the governing equations according to the pressure jump at the 
cell interface [15]. Note that any discontinuity in the reconstructed data means the cell resolution is 
not enough to resolve the flow structure. Here, the evaluation of the critical Rayleigh numbers for the 
two-dimensional compressible convection case is a good indication of the accuracy and robustness of 
the scheme. 

In this section, we adopt the cases calculated by Grahma [6] and Chan [3] for the purpose of 
comparison. Consider an ideal gas confined in a box with < xi < /, < X2 < d. The free surface 
and isothermal boundary conditions are imposed at the top and bottom boundaries, 

U2 = 0, ^=0, T = Tb,Tt, at X2 = 0,d, 

dX2 

and the stress-free insulated boundary conditions are imposed in the horizontal direction, 

^ du2 dT ^ 

ui = 0, - — = 0, - — = 0, at xi = 0, /. 

dxi dxi 

The initial condition is the static state with a small perturbation in the velocity field, 

T = {l + Z{d-X2)/d)Tt, 
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p = {T/Ttr+'pt, 

where Z = {Tj, — Tt)/Tt is the normahzed parameter, and m is the polytropic gas index. There is a 
hydrostatic solution of the Navier-Stokes equations with a uniform constant gravitational acceleration, 

d 

With the above T, p, and p distributions, the initial vertical velocity is perturbed by 

ui = 0, 

. 2ttxi . 7rx2 
U2 = «o sm — ; — sin — — , 
/ d 

where uq is a very small constant. Actually, the initial perturbation for the vertical velocity may 
be randomly generated with the satisfaction of the boundary condition. The above form leads to a 
faster startup of the convective rolls. All the quantities, pt, pt, Tt, d and I, are normalized. In all 
computations presented here, 7 = 5/3, Pr = 1 and m = 1.4 are used. The Rayleigh number is defined 
as 



1 — (7 — l)m 



(m + 1) 



7 

The standard solutions of the two dimensional laminar convection are shown in Figure 6. A discussion 
of these solutions can be found in [4]. 

For each Z, the calculation is repeated with two slightly supercritical Rayleigh numbers, Rai and 
Ra2- Then, the critical Rayleigh number is determined by fitting the curve max (^2) = A exp {B(Ra — Ruc 
by a least square approximation, where U2 is the vertical velocity, Rac is the critical Rayleigh number, 
and A and B are the constants to be determined. The results given in Table 1 are obtained from 
the current BGK scheme, where a continuous interpolation (4) is used in the reconstruction stage. 
From Table 1, we can see that the current BGK scheme is an accurate method for a viscous flow 
under a gravitational field. The difference between the critical Rayleigh numbers calculated from the 
gas-kinetic BGK scheme and those from the linear analysis is less than 3.5% for Z < 10 and < 1% for 
Z < 4. Also, the Courant number used for each calculation is included in Table 1. When a nonlinear 
limiter, Eq.(5) and Eq.(6), is used in the reconstruction of the initial data, discontinuities at the cell 
interfaces will be introduced. This will generate extra kinematic dissipation. Certainly, the interface 
discontinuity will become smaller for higher-order interpolation in the smooth region, as will the nu- 
merical dissipation. But, instead of using higher-order interpolation, the van Leer limiter is always 
used in this paper because of its robustness and simplicity. To investigate the effect of the limiter 
on the Rayleigh number, the gas-kinetic BGK scheme with the van Leer limiter is used for the same 
calculation. The numerical computation shows that, for the weakly compressible flow, the critical 
Rayleigh number becomes slightly larger {Ra = 553.67 for Z = 0.5), while for highly compressible 
flow, the situation is complex. The critical Rayleigh number becomes a little bit smaller {Ra = 289.51 
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) for Z = 2 and larger {Ra = 167.5) for Z = A. For Z = 10, the instability appears near the upper 
surface. This can be explained by the use of a limiter, which continuously perturbs the boundary layer 
through the introduction of numerical discontinuities. The limiter generates a jump in the vertical 
velocity field. A direct remedy for this is to use the high-resolution interpolation, or modify the initial 
dissipative terms in /o to make these terms continuous across the cell interface [10]. However, this 
error due to the limiter could be ignored for small and modest Z. In conclusion, a nonlinear limiter 
will introduce numerical dissipation, but its effect on the calculation of the critical Rayleigh number 
is tolerable. Furthermore, for Pr <^ 1 cases, the effect from the nonlinear limiter becomes negligible, 
where the heat conduction dominates and smoothes the flow field. 

3.4 Thermodynamic properties of the laminar convection in 3D 

In the numerical flux, the derivatives of the flow distribution in the tangential direction of a cell inter- 
face are taken into account in the current scheme. The current scheme is essentially a multidimensional 
method, where the physical effect from both gradients in the normal and tangential directions is taken 
into account. In order to verify the necessity of this kind of discretization, tests of compressible 
convection in a three-dimensional box will be presented in this subsection. 

The two-dimensional model in the previous subsection is extended to three dimensions. The aspect 
ratio of the rectangular box is 1:0.1:1. The insulated solid well boundary conditions are applied in 
the horizontal directions (xi and X2). Initially, the velocity field is perturbed in all directions: Xi, X2 
and 2:3. The Courant number takes a value of 0.3 and the grid size is 50 x 5 x 50. Unless otherwise 
specified, all the parameters in this subsection have the same value as those in the previous subsection. 

Among the solutions, the maximum Mach number and Nusselt number are used to measure the 
thermodynamical property and convective heat transfer. The maximum Mach number is defined as 



Mc = max 



which sensitively depends on the viscosity. 

The Nusselt number is used to compare the flux carried by the heat conduction and by convection. 
Following Graham [6], it is defined as. 



F — F 



where F^, = —/S.(j)k/cp, = K{Ti,—Tt)/d are the heat flux due to the adiabatic gradient and conduction, 
respectively. Let 

Ft = pCpTus + K— 

be the total heat flux, where the over-bar means the average over the whole computational domain. 
Numerical experiments show that M cannot be correctly obtained with a low order scheme. We also 
introduce the relative Rayleigh number: 

where R is defined by equation (20) and Rc is the critical Rayleigh number from the linear analysis. 
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A series of finite-amplitude steady solutions is calculated with different Rayleigh numbers. The 
depth parameter is fixed to a unit. The results of one-cell circulation are given in Table 2 and Table 3. 
The comparison with Graham's results is shown in Figure 7. 

Prom plot (a) in Figure 7 we can see that the agreement is good when the central interpolation is 
used for the initial reconstruction (the average deviation a < 9%). The agreement is better for larger 
Rayleigh numbers, for example, (j{R* = 100) < 5%. This is expected since a BGK scheme with a 
central interpolation is similar to the Lax-Wendroff scheme, which was used by Graham [6]. As pointed 
out in the previous section, the numerical dissipation will be introduced when the van Leer limiter is 
adopted. The average deviation of the results with different reconstructions is less than 2%. This kind 
of discrepancy increases with the Rayleigh number {cr{R* = 3) = 1.4%) and a{R* = 100) = 3.4%). 

The trend of the discrepancy in plot (b) of Figure 7 is similar to that in plot (a) , except that the 
best fitting results come with the van Leer limiter. The average deviations from Graham's results 
are approximately less than 3% and 5% using the van Leer limiter and the central interpolation, 
respectively. The Nusselt numbers obtained from the BGK scheme are systematically larger than 
Graham's results. 

The effects of including tangential derivatives, i.e., the so-called multi-dimensionality, into the 
numerical fluxes have been checked. The results obtained from a dimensional-splitting BGK method, 

c7"+i = if,3c7,"+^ + Q(?7r+i), 

arc also listed in Tables 2 and 3, where is the one-dimensional operator and Q is the source term. 
In the computation of the numerical fluxes, the tangential derivatives are excluded. In the directional 
splitting method, it is clear that additional numerical dissipation is added to the scheme because the 
maximum Mach numbers in the third row of Table 2 are systematically less than those in the second 
row. The situation is complicated for the Nusselt number. When a dimensional splitting method 
is used, the convection becomes slightly turbulent for R* > 70. The numbers with a superscript 
* represent temporally averaged values. The standard deviation for Mc{R* = 100), N{R* = 70) 
and N{R* = 100) are 0.001, 0.003 and 0.226, respectively. When the flow becomes turbulent, the 
Nusselt numbers from the multidimensional scheme are a bit smaller than those from the dimensional 
splitting scheme which indicates that the turbulence has a more efficient way to transfer the heat 
than the laminar flow. This test clearly shows that the multidimensional scheme is more stable and 
accurate than the dimensional splitting scheme. The reason lies in the fact that to keep both normal 
and tangential derivatives in the evaluation of a flux function in a circulation motion is close to the 
physical reality. 

3.5 Rayleigh- Taylor instability 

This test case is similar to LeVeque and Bale's Rayleigh-Taylor instability case [8], which was used 
to validate their quasi-steady wave propagation scheme. Consider an isothermal equilibrium ideal gas 
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(7 = 1.4) in a 2D polar coordinate, 



Poir) 



_ -a{r+ro) 



a 

-V(j){r) = -1.5, 
no = 0, 



where a = 2.68, ro = 0.258 for r < ri and a = 5.53, ro = —0.308 for r > ri with 

n = 0.6(1 + 0.02 cos (20^)) 

for density and 

n = 0.62324965 

for pressure. In this case, two different isothermal equilibrium states intersect around r = 0.6 with 
density and pressure perturbations. 

Due to the geometric symmetry in this problem, only the solution in the first quadrant is computed. 
The computational domain, [0, 1] x [0, 1], is covered by 120 x 120 uniform grid points. The CFL number 
takes a value 0.6. On the axis of symmetry, a rotated periodic boundary condition is used. At the other 
boundaries, the values in the ghost cells are fixed as their initial values in order to keep the isothermal 
equilibrium solution there. Figure 8 presents the time evolutions of the density distributions at times 
t = 0.0, 0.8, 1.4, and 2.0. Figure 9 shows a scatter plot of the density as a function of the radius. Note 
that, away from the location (r ~ ro), the current method retains the hydrostatic equilibrium solution 
and radial symmetry very well. 

4 Conclusions 

The gas-kinetic BGK scheme for the Navier-Stokes equations is extended to include external gravi- 
tational forces. Our study allows us to draw the following conclusions. (1) The gravitational force 
effect should be included in the evaluation of the numerical fluxes across a cell interface. It has the 
first-order effect in time [16], which cannot be neglected in a second-order scheme. In a gas kinetic 

scheme, the external force modifies the time evolution of both the initial gas distribution function and 
the local equilibrium state through the particle acceleration. Without such an implementation, it is 
impossible for the system to settle down to a highly accurate hydrostatic equilibrium state. (2) Nu- 
merical heating is mainly caused by the inconsistent treatment of the gravitational source terms inside 
each cell in the energy source term scheme. After a long time evolution, the numerical error will accu- 
mulate and dramatically affect the final solution of an isolated or insulated gravitational system. For 
a system that exchanges energy with its surrounding environment or that has a short evolution time, 
such as the star explosion case, this kind of numerical heating may not appear to be significant. (3) 
The numerical simulation of the onset of two-dimensional laminar convection shows that the current 
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BGK scheme can accurately determine the critical Rayleigh numbers. However, the nonlinear hmiter 
is usuahy harmful to the smooth flow solution. High-order interpolation may alleviate this problem. 
(4) The results of finite amphtude laminar convection with a central interpolation agree weh with 
Graham's results. The van Leer hmiter introduces a bit of additional dissipation. By including the 
tangential gradient terms, i.e., the so-called multidimensional approach, the accuracy of the scheme 
is improved. The kinetic scheme presented in this paper is not only limited to hydrodynamics under 
gravitational fields; it can be applied to other systems in which the external force term can be an 
electric field force, such as the Boltzmann-typc equations for semiconductor devices. 

Developing a well-balanced scheme for gas dynamic equations under gravitational fields is still an 
open problem. It is much more difficult than developing a scheme for the shallow water equations. 
We hope that this paper will help to shift attention from shallow water equations to gas dynamic 
equations. At the current stage, it seems that no currently available scheme can maintain the velocity 
of the hydrostatic equilibrium state up to the machine zero, such as on the order of 10~^^. At the same 
time, the scheme has a shock capturing ability. Even with consideration of all the techniques presented 
in this paper, the current method has only second-order accuracy in maintaining the hydrostatic 
solution. Much work is left for the future. 
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Appendix: Solution of Matrix Equation b = Ma in the three-dimensional 
Case 

In the gas-kinetic scheme, to obtain the derivatives of a Maxwellian, the solution of the following 
equations need to be evaluated. 



(h 
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where b and M are known. The matrix M is defined as 

M = - / ^0 i^f^dE, 
PJ 

where is the Maxwellian. 

The details of the matrix M are given by 
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where 
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-iul + ul + ui + iK + 3)/2X), 



-{ul + ui{ul + ul) + {K + 5)ui/2X), 
-{ul + U2{ul + ul) + {K + 5)U2/2X), 
-{ul + U3{ul + ul) + {K + 5)u3/2\), 



and 



B5 = -i{{ui +ui + uD' + {K + 5){ui + ui + ui)/X + {K + 3){K + 5)/4A' 



Defining 



R5 = 265 

i?4 = b4-U3bi, 

R3 = 63 -"2^1, 

R2 = h — uibi, 



2 , 2 I 2,-^ + 3 
U1+U2+U2 + 



2A 



^1, 



tiie solutions of Eq.(21) are 



4A2 



05 
03 



K + 3 
2XR4 — 1*305, 



(i?5 - 2uiR2 - 2U2R3 - 2U3R4), 



2XR3 - 142(15, 

02 = 2Ai?2 - uia^, 



and 
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Table 1: Critical Rayleigh number. 



z 


0.5 


1 


2 


4 


10 




545.03 


422.79 


288.51 


164.88 


58.17 




550.02 


424.84 


289.82 


167.24 


59.96 


5" 


0.4 


0.2 


0.1 


0.07 


0.02 



Linear analysis results from Graham [6]. 
Results from current gas-kinetic BGK scheme. 
Courant number. 



Table 2: Maximum Mach number (Mc) vs. relative Rayleigh Number (R*). 



R* 


3 


7 


10 


30 


70 


100 


Mc'' 


0.072 


0.097 


0.107 


0.136 


0.161 


0.172 




0.071 


0.096 


0.106 


0.134 


0.156 


0.166 




0.067 


0.095 


0.104 


0.131 


0.151 


0.159* 



Current BGK scheme with central interpolation. 
Current BGK scheme with the van Leer limiter interpolation. 
Dimensional-splitting BGK scheme with the van Leer limiter. 
Temporally averaged value. 
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Table 3: Nusselt number (AA) vs. relative Rayleigh Number 



R* 


3 


7 


10 


30 


70 


100 




2.575 


3.729 


4.334 


6.938 


10.000 


11.676 




2.533 


3.702 


4.299 


6.831 


9.675 


11.141 




2.236 


3.393 


3.989 


6.777 


9.849* 


11.252* 



Current BGK scheme with central interpolation. 
Current BGK scheme with the van Leer limiter. 
Dimensional-splitting BGK scheme with the van Leer limiter. 
Temporally averaged values. 
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Figure 1: Numerical heating. The time evolution of the total internal energy of an isolated system. 
Both solutions are computed by the BGK scheme but with different source term treatment. The total 
internal energy is scaled by the initial value. The dashed line is the solution from the BGK scheme 
without including the source term effect in the flux function and the solid line is the solution with the 
source term effect included in the flux. 
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Figure 2: Solutions of the BGK scheme with and without the source term effect in the flux function. 
The BGK scheme with the van Leer limiter is used in these calculations. 100 grid points are used in 
the whole domain. The output time is t = 0.25. The upper and lower figures have the perturbations; 
upper: 77 = 0.01, lower: rj = 0.001. The dotted lines show the initial pressure perturbation from a 
hydrostatic solution, the pluses represent the results without the source term effect in the flux, and 
the solid lines are the results with the source term effect. The diamonds are the reference solutions 
from [8]. 
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Figure 4: Results of the ECT scheme with the external force term in the flux (diamonds), and the ECT 
scheme without including external force term in the flux (solid lines). These results were obtained 
after 500, 000 time steps. 
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Figure 6: Solutions of two-dimensional laminar convection obtained by the BGK scheme. The upper- 
left, upper-right, lower-left and lower-right plots show the velocity vector field, the density distribution, 
the pressure distribution, and temperature distribution, respectively. The corresponding parameters 
are Z = 10 and Ra = 61. 
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Figure 7: Maximum Mach number (a) and Nusselt number (b) vs. relative Rayleigh numbers. The 
solid lines are the results from Graham's paper [6] and the diamonds and pluses represent the results 
computed by the BGK schemes with central interpolation and the van Leer limiter reconstruction. 
The triangles show the results from the dimensional-splitting BGK scheme with the van Leer limiter. 
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Figure 8: A Rayleigh- Taylor instability with gravitational field directed radially inward. Density 
contours at four different times {t = 0.0,0.8,1.4, and 2.0) are shown in the four quadrants, starting 
with the initial data in the upper right corner and progressing clockwise. 
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density scatter at t = density scotter ot t=0.8 




Figure 9: Scatter plots of the density in the cell vs. the distance of the cell center from the origin. 
The density is shown at four different times {t = 0.0,0.8, 1.4, and 2.0). The single line away from ro 
indicates the capacity of the current scheme to keep the hydrostatic equilibrium solution. 
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